% By Martin M. Andreasen, Auguest 2009
% This function computes the conditional moments up to k_period into the
% future for the variables selected by "y_select" in the g-function. 
% This implementation allows for non-zero third moments in the shock
% distributions.
% Note that y_select must have the form of 1 x dimy, where dimy is the
% number of variables in the g-functions we want to compute conditional
% moments for. The steady state value of y_select should be stored in y_ss
%
% In terms of the notation we solve: p(x,sig) = E_t[r(x_t+1,sig)]
% The computations are derived by using the Perturbation On Perturbation (POP) method.
% This function is for a log-transformation
%
function [px,pxx,pss,pxxx,pssx,psss] = CondMoments_Mom3_log(gx,gxx,gss,gxxx,gssx,gsss,hx,hxx,hss,hxxx,hssx,hsss,eta,vectorMom3,y_select,k_period,order_app);

% Allocating memory
dimy    = size(y_select,2);
nx      = size(hx,1);
ne      = size(eta,2);
px      = zeros(dimy,k_period,nx);
pxx     = zeros(dimy,k_period,nx,nx);
pss     = zeros(dimy,k_period);
pxxx    = zeros(dimy,k_period,nx,nx,nx);    
pssx    = zeros(dimy,k_period,nx);
psss    = zeros(dimy,k_period);

for i=1:dimy
    rx   = reshape(gx(y_select(1,i),:),1,nx);
    rxx  = reshape(gxx(y_select(1,i),:,:),nx,nx);
    rss  = gss(y_select(1,i),1);
    rxxx = reshape(gxxx(y_select(1,i),:,:,:),nx,nx,nx);
    rssx = reshape(gssx(y_select(1,i),:),1,nx);
    rsss = gsss(y_select(1,i),1);
   
    for j=1:k_period
       % ************** The first order effects *****************
       px(i,j,:) = rx(1,:)*hx;
       
       
       %*************** The second order effects ***************
       if order_app > 1
           tmp = zeros(nx,nx);
           for gama1=1:nx
               tmp = tmp + rx(1,gama1)*squeeze(hxx(gama1,:,:));
           end
           pxx(i,j,:,:) = hx'*rxx*hx + tmp;

           pss(i,j) = rx(1,:)*eta*eta'*rx(1,:)' + trace(eta'*rxx*eta) + rx(1,:)*hss + rss;
       end
       
       % ***************** The third order effects *****************
       if order_app > 2
           for alfa1=1:nx
               for alfa2=1:nx
                   tmp = zeros(1,nx);
                   for gama3=1:nx
                       tmp = tmp + hx(:,alfa1)'*reshape(rxxx(:,:,gama3),nx,nx)*hx(:,alfa2)*hx(gama3,:);
                   end
                   pxxx(i,j,alfa1,alfa2,:) = tmp  ...
                       + hx(:,alfa1)'*rxx*reshape(hxx(:,alfa2,:),nx,nx) ...
                       + hx(:,alfa2)'*rxx*reshape(hxx(:,alfa1,:),nx,nx)...
                       + reshape(hxx(:,alfa1,alfa2),1,nx)*rxx*hx(:,:)...
                       + rx(1,:)*reshape(hxxx(:,alfa1,alfa2,:),nx,nx);
               end
           end

           pssx(i,j,:) = 2*rx(1,:)*eta*eta'*rxx*hx+hss'*rxx*hx+rx(1,:)*hssx+rssx*hx;
           for gama3=1:nx
               pssx(i,j,:) = reshape(pssx(i,j,:),1,nx) + trace(eta'*squeeze(rxxx(:,:,gama3))*eta)*hx(gama3,:);
           end

           tmp = 0;
           for phi1=1:ne
               tmp = tmp + rx(1,:)*eta(:,phi1)*rx(1,:)*eta(:,phi1)*rx(1,:)*eta(:,phi1)*vectorMom3(phi1)...
                   + 3*eta(:,phi1)'*rxx*eta(:,phi1)*rx(1,:)*eta(:,phi1)*vectorMom3(phi1);
           end
           for gama1=1:nx
               for phi1=1:ne
                   tmp = tmp + eta(:,phi1)'*reshape(rxxx(gama1,:,:),nx,nx)*eta(:,phi1)*eta(gama1,phi1)*vectorMom3(phi1);
               end
           end
           tmp = tmp + rx(1,:)*hsss + rsss;
           psss(i,j) = tmp;
       end

       % Updating rx
       rx   = reshape(px(i,j,:),1,nx);
       rxx  = reshape(pxx(i,j,:,:),nx,nx);
       rss  = pss(i,j);
       rxxx = reshape(pxxx(i,j,:,:,:),nx,nx,nx);
       rssx = reshape(pssx(i,j,:),1,nx);
       rsss = psss(i,j);
    end
end

% We reexpress the output if dimy = 1
if dimy == 1
    px      = squeeze(px(dimy,1:k_period,1:nx));
    pxx     = squeeze(pxx(dimy,1:k_period,1:nx,1:nx));
    pss     = reshape(pss(dimy,1:k_period),k_period,1);
    pxxx    = squeeze(pxxx(dimy,1:k_period,1:nx,1:nx,1:nx));    
    pssx    = squeeze(pssx(dimy,1:k_period,1:nx));
    psss    = reshape(psss(dimy,1:k_period),k_period,1);
end